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Abstract. Aiming for simplicity of explicit equations and at the same time 
controllable accuracy of the theory we present results for all thermodynamic 
quantities and correlation functions for the weakly interacting Bose gas at 
short-to-intermediate distances obtained within an improved version of Bcliaev's 
diagrammatic technique. With a small symmetry breaking term Beliacv's 
diagrammatic technique becomes regular in the infrared limit. Up to higher-order 
terms (for which we present order-of-magnitude estimates), the partition function 
and entropy of the system formally correspond to those of a non-interacting 
bosonic (pseudo-)Hamiltonian with a temperature dependent Bogoliubov-type 
dispersion relation. Away from the fluctuation region, this approach provides 
the most accurate — in fact, the best possible within the Bogoliubov-type pseudo- 
Hamiltonian framework — description of the system with controlled accuracy. 
It produces accurate answers for the off-diagonal correlation functions up to 
distances where the behaviour of correlators is controlled by generic hydrodynamic 
relations, and thus can be accurately extrapolated to arbitrarily large scales. In 
the fluctuation region, the non-perturbative contributions are given by universal 
(for all weakly interacting U(l) systems) constants and scaling functions, which 
can be obtained separately — by simulating classical U(l) models — and then used 
to extend the description of the weakly interacting Bose gas to the fluctuation 
region. The theory works in all spatial dimensions and we explicitly check its 
validity against first-principle Monte Carlo simulations for various thermodynamic 
properties and the single-particle density matrix. 
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1. Introduction 

For nearly half a century the theory of the weakly interacting Bose gases (WIBG) 
remained in the realm of purely theoretical investigations [TJ [2], SI [3l [6J [T] (for a recent 
review, see [5]) providing insight into the nature of superfluid states of matter but not 
directly relating to existing experimental systems. The situation changed with the 
realization of Bose Einstein condensation in cold atomic gases [9l [KB [TT] . The typical 
values of the interaction parameter for alkali atoms are very small na 3 ~ 10~ 6 , where 
n is the number density of the gas and a is the s-wave scattering length. Since 1995 
both the theory of WIBG and experiments have progressed with a close relationship 
between the two [12] . 

Until recently, existing mean-field and variational treatments such as the 
Bogoliubov zero-temperature approximation and the finite-temperature quasi- 
condensate theory [13j were capable of describing the data within the relatively large 
experimental uncertainties. However, the improvements in detection techniques, the 
studies of optical lattice systems in which the effective gas parameter can be made 
arbitrarily large, and the need for reliable thermometry, e.g., through precise entropy 
matching, all indicate the need for a more accurate and controllable theory. It is 
highly desirable, similarly to the well-known corrections to the ground state energy 
[T4] and condensate density [I] , to have a theory which accounts for leading corrections 
to all thermodynamic properties at finite temperature, works in any dimension, and 
provides estimates for omitted higher order terms. 

The accuracy of the description plays an important role in this paper. We 
systematically estimate the errors of all approximations, which is crucial for comparing 
seemingly different schemes proposed for Bose gases. Indeed, alternative theories can 
be equivalent to each other within the level of their accuracy, but a definitive conclusion 
can only be drawn on the basis of analyzing their systematic errors. For example, 
within the same accuracy gap-less and gap-full finite-temperature schemes may often 
be regarded as equivalent without any preference for one of them, as long as the gap 
is on the order of or smaller than the terms in the chemical potential which the theory 
ignores. Also, in the fluctuation region all perturbative schemes are equally inaccurate 
in terms of their treatment of the order parameter, and it makes no sense whatsoever 
to distinguish between them by the type of the transition they predict. 

In this work, we provide a rigorous framework for obtaining a consistent 
description of all thermodynamic finite-T properties of the WIBG in one, two, and 
three dimensions. It is based on Bcliaev's regularized diagrammatic technique [2] 
which allows us to calculate all relevant thermodynamic functions of the system. The 
diagrammatic technique also leads to an accurate description of correlation functions 
up to sufficiently large distances where Popov's hydrodynamic description takes over 
[5]. The key results of our calculations are summarized in section[2j which is accessible 
to the general reader. 

It turns out that all final results for thermodynamic functions perfectly fit into the 
pscudo-Hamiltonian picture, i.e. we prove that the same answers would be obtained 
if they were calculated in the standard way with Bogoliubov's non-interacting quasi- 
particle picture (with a self-consistently defined temperature dependent spectrum). 

In addition, we find an explicit expression for the pressure and give a more 
accurate expression for the chemical potential. We obtain the equations of state for 
all basic thermodynamic quantities in an integral parametric form using temperature, 
effective chemical potential, and the interaction strength as parameters. These 
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Figure 1. Sketch of the density matrix dependence on distance in the WIBG at 
temperatures below the fluctuation region. After an initial drop the density matrix 
levels-off at the healing length and undergoes slow decay up to a distance where 
the supcrfiuid hydrodynamic approach becomes applicable. The U(l) symmetry 
breaking terms make the density matrix to converge to a finite value (in any 
dimension) at the hydrodynamic length scales. 

results include sub-leading corrections [J The analysis of systematic errors shows 
that the accuracy of our results cannot be further improved within the paradigm of 
non-interacting quasiparticlcs with a (temperature-dependent) Bogoliubov spectrum. 
Moreover, there exists a two-parametric continuum of possibilities for recasting the 
form of the final answers with the same accuracy due to the freedom of modifying the 
definitions of approximate notions of "effective interaction" and "effective chemical 
potential" . Using this freedom, one can trade the higher order terms in the explicit 
integral expressions for the thermodynamic quantities for those of the effective 
interaction vertex. Our analysis shows that the question of the energy and momentum 
dependence of the interaction vertex is more a matter of taste rather than accuracy. 
We show that even in 2D — where, in view of the vanishing s-scattering amplitude, a 
low-energy cutoff is unavoidable and one could expect that the momentum-dependent 
effective-interaction vertex is crucially important in any theory taking care of sub- 
leading corrections — the effects associated with the energy and momentum dependence 
of the effective interaction can be naturally accounted for on equal footing with all 
the other effects of the same order. Within the systematic error bars of our pseudo- 
Hamiltonian approach, the effective interaction of the 2D WIBG can be safely treated 
as an energy and momentum independent constant. 

Such analytical approaches are not supposed to work in the vicinity of the critical 
point T c . This so-called fluctuation region is characterized by scaling functions that 
are universal for all U(l) weakly-interacting systems. It has already been studied 
numerically with high accuracy |15[ 1161 117j , so that a combination of analytical and 
numerical data provides an accurate description at all temperatures. 

\ The only exceptions are the quantities, like entropy and heat capacity, that vanish as T — * 0. At low 
enough temperatures, these quantities lose their ideal-gas leading terms, getting nothing instead — as 
opposed to energy, chemical potential, and pressure, having straight-forward mean-field contributions 
that start to lead at low temperatures. 
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Despite its success in treating the weakly interacting Bose gas, Bcliaev's original 
technique features a subtlety that has seriously questioned its reliability: the normal 
and anomalous self energies — the fundamental quantities of the theory — have been 
demonstrated to be essentially momentum-dependent, with the anomalous self-energy 
vanishing in the limit of zero momentum and frequency [TSl [TH [20] . On the other hand 
the first-order diagrams in the theory of the weakly interacting gas with a contact 
interaction yield momentum independent self energies. In the literature, one can find 
a number of recipes of how to deal with this inherent infrared problem of Beliaev's 
technique. In Popov's approach [5], the diagrammatic technique applies only to the 
higher-momentum part of the bosonic field, while the lower-momentum part is treated 
separately within the hydrodynamic representation. Similar results can be achieved by 
working with a finite-size system and taking advantage of the bi-modality of correlation 
properties of weakly interacting superfluid bosonic systems (away from the fluctuation 
region) [2T]. The bi-modality allows one to select such a system size that in terms of 
local thermodynamic properties the system is already well in the macroscopic limit 
while the effect of the infra-red renormalization of the self energies is still negligibly 
small. Another option of a controllable microscopic description of weakly interacting 
Bose gas is the renormalization-group treatment [22\ 123] . Similar to the finite-size 
regularization, our approach is to exploit the bi-modality of correlation properties 
of the WIBG, see figure [TJ by cutting off infra-red divergencies using small U(l) 
symmetry breaking terms. The amplitude of the symmetry-breaking terms is such that 
their effect on the local thermodynamic properties of the system is negligibly small 
while all infrared singularities in Beliaev's technique are gone. Clearly, the solution to 
the infrared problem comes at the expense of suppressing long-range fluctuations of 
the phase of the order parameter (with opening of a gap in the spectrum at small 
momenta), and thus distorting the long-range behaviour of correlation functions. 
Nevertheless, this distortion is essentially irrelevant since it takes place at large 
distances where the (generic to all superfluids) behaviour of correlators is governed 
by hydrodynamic fluctuations of the phase of the order parameter, with a simple 
and transparent theoretical description |24j . Moreover, in all final answers one can 
explicitly take the limit of vanishing symmetry breaking terms. 

In terms of the final results, we confine our analysis to the most typical (and 
most difficult in two and three dimensions) case of the dilute gas, when the size of the 
potential Rq is much smaller than the distance between the particles: 

i?o«"~ 1/d . (1.1) 

In principle, the theory works also at Rq > n _1 / d , and becomes even simpler, since in 
this case the weakness of interaction literally implies Born type of the potential (while 
in the opposite case the amplitude of potential can be arbitrarily large in two and 
three dimensions). However, only with the condition (|f .![) the final answers become 
insensitive to the details of the interaction potential. The second condition which is 
assumed in our final results is 

At » Ro, (1.2) 

where At is the de Broglic wavelength. Given the inequality (|1.1|) , the condition (| 1 . 2[) 
is satisfied automatically as long as one is interested in essentially quantum properties 
of the system, implying X^n > 1. An extension of the theory to the Boltzmann high- 
temperature regime At ^ Ro-, where answers become sensitive to the particular form 
of the interaction potential, goes beyond the scope of this paper. 
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The paper is organized as follows. In section[2]we first summarize the key results 
of our approach for thermodynamic quantities and correlation functions in a form 
easily accessible to the general reader. The next chapters systematically introduce the 
technique and derive the results. In section [3] we re-derive the diagrammatic technique 
for the phase with the broken U(l) symmetry by introducing explicit (arbitrary small) 
symmetry breaking terms in the Hamiltonian. We keep the discussion as general 
as possible to cover inhomogeneous, e.g. trapped, systems and states with non-zero 
average current. In section [3] we also derive the central expressions for the proper self 
energies and chemical potential. In section 3] we consider thermodynamics (including 
the superfluid density) of a homogeneous system in the superfluid phase. In particular, 
we introduce the notion of the renormalized chemical potential fi and show that 
fj, and temperature form a natural pair of thermodynamic parameters, in terms of 
which all the other quantities are obtained from the single-particle momentum-space 
integrals. In section Owe discuss, at the order-of-magnitude level, the structure of 
higher-order corrections to the spectrum of elementary excitations and thermodynamic 
quantities. This analysis, in particular, leads to the conclusion that within the pscudo- 
Hamiltonian approximation and Bogoliubov ansatz for the spectrum of elementary 
excitations our results cannot be further improved. In section [6] we show that the 
answers for the asymptotic long-range behaviour of off-diagonal correlation functions 
are readily obtained on top of Bcliaev's diagrammatic treatment for medium-range 
correlations, by employing the generic hydrodynamic description of a superfluid. The 
short section [7] deals with the normal region. In section[8]we show that in both normal 
and superfluid regimes (but away from the fluctuation region around the critical 
point) the thermodynamic relations can be associated with a pseudo-Hamiltonian 
in the following sense. The partition function and the entropy of the system turn out 
to correspond to those of a non-interacting bosonic Hamiltonian with temperature- 
dependent parameters and a temperature-dependent global energy shift. In section 
[9] we render the results for the fluctuation region, which previously were obtained 
in 2D and 3D, but not in ID. In section [10] we finally compare our solutions with 
first-principle Monte Carlo data for WIBG both in continuous space and on a lattice. 

2. Key results 

Here we summarize the most important relations which we derive in the main part 
of the paper, assuming that the conditions (|1.1|) and (|1.2[) are met. The results of 
this section work for both continuous-space and lattice systems; in the latter case, the 
mass to is understood as the effective mass of low-energy motion. The equilibrium 
state of the system is conveniently parametrized by two independent variables, the 
temperature T and the effective chemical potential fi. The latter is always negative, 
so that fi = —\p\. 

2.1. Coupling constant 

The interaction is characterized by an (effective) coupling constant U. In one 
dimension U is the zero- momentum Fourier component of the bare potential U(r). In 
two and three dimensions, U is naturally expressed in terms of the s-wave scattering 
length a; the latter being defined as the radius of hard disk/sphere potential 




r < a , 
r > a , 



(2.1) 
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with low-energy scattering properties identical to those of the potential U(r). A 
representational freedom allows one to use slightly different expressions for U , without 
sacrificing the order of accuracy. Our particular choice is 

AlTfl 

U = (d = 3), (2.2) 

TO 

in three dimensions, and 

U= , (9/ t /m , o ( d = 2 )> ( 2 - 3 ) 
m(2/a z meo) — 2j 

in two dimensions, where 7 = 0.5772 ... is Euler's constant. A particular expression for 
U comes from fixing the specific form of an auxiliary function II(fc). The latter enters 
thermodynamic relations in such a way that the values of those remain insensitive — up 
to irrelevant higher-order corrections — to slightly changing U by adopting a slightly 
different II. Equations (|2.2[) and (|2.3[) correspond to 

nW=-^ (<*=3), (2.4) 

respectively, where e(k) is the particle dispersion relation. In one dimension II = 
0. Within the representational freedom, different choices of the (II, J7)-pair are 
straightforwardly connected with each other by (|3.49|) - (|3.50p . 

In the 2D case, the coupling constant U, and, correspondingly, function Il(fc), 
slowly evolve with changing effective chemical potential and temperature: 

..«,(M0-{ If I*' f T lf- (2.6) 



It is well within the representational freedom to multiply the r.h.s. of (|2.6|) by a 
constant of order unity, provided the same expression for eo is used in both (|2.3p and 
(|2.5[) . The only reason for introducing the order-unity factor e/2 is to cast the zero- 
temperature equations of state in the form (|4.46[) - (|4.47[) adopted previously by some 
authors. 

In the case when the problem of finding the scattering length a is not 
straightforward, one can resort directly to the integral equation (|3.41[) defining U 
in terms of U and II. Numeric solution of this equation can be obtained within bold 
diagrammatic Monte Carlo technique |25j . which gives the scattering length as a by- 
product of the calculation. 

2.2. Thermodynamics in the super fluid region 

Thermodynamic quantities are obtained in parametric form, as functions of the pair 
of parameters jl and T. The main relation is for the number density: 

n(fi,T) = \j2\/U + n', (2.7) 



e(k)-E(k) , , , N e(k) 



(2.8) 



where 



E(k) = y/e(k)[e(k) + 2\ii\] (2.9) 
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is the Bogoliubov-type quasiparticlc dispersion relation and 

-l 



E 



,E/T 



1 



(2.10) 

is the Bose-Einstein occupation number. In equation (|2.8[) and below throughout 
the text where we sum over momenta, we set the system volume equal to unity thus 
omitting a trivial normalization factor. 

The genuine chemical potential is then given by 

fj, = nU — |/2| . 

The pressure and entropy density are obtained as 

p = p 2 /2U + 2n'\p\ + (n) 2 U 



-E(k)/T 



(2.11) 



}, (2-12) 



and 



E 



^N E In I 



(2.13) 



from which the energy density follows from the general relation e = fxn + Ts — p. 

At zero temperature the integrals can be evaluated analytically, as discussed in 
subsection 14.41 

In continuous space, where e(fc) = k 2 /2m, the superfluid density is given 
compactly by Landau's formula 

1 ^ k 2 dN , AS 

n s = n+- V — j—. 2.14 

k 

Due to the condition (jl.ljl the answer for a lattice differs only by a global factor equal 
to the ratio between the bare and effective masses. 



2.3. Single-particle density matrix. Quasi- condensate 



In the superfluid region, the single-particle density matrix can be parameterized as 



p(r) = p(r)e- A «, 



with 



p(r) 



k 



1 - e* 



r 

E 2 



E-e-|A| 



+ {E + \p,\)N E 



A(r) = M£(l- 

T7 * v 



,ikr 



E-e 
~2E 2 ~ 



[1 + 2N E ] 



(2.15) 



(2.16) 



(2.17) 



By definition, p(oo) is the condensate density. The function p(oo) is finite in any 
dimensions. In ID and at finite temperature in 2D, the function A(r) diverges at 
r — > oo, consistently with the general fact of absence of condensate in those systems. 
Nevertheless, at distances at which the function p(r) saturates to its asymptotic value 
p(oo), the function A(r) is still much smaller than unity. This characteristic feature 
of the weakly interacting system allows one to speak of the quasi-condensate with the 
density given by p(oo). Up to the distances at which A becomes ~ 1, the correlation 
properties of condensed and quasi-condensed systems are essentially the same. Details 
of the general approach to long-range off-diagonal many-particle correlation functions 
arc described in section [6] 



Beliaev technique for a weakly interacting Bose gas 



8 



2-4- Thermodynamics in the normal region 

The density n = n(fl, T) in the normal region is given by 

n = E [ e ' (fe)/T - l ] _1 > ( 2 - 18 ) 

k 

with 

e(fc) = e(fc) + |£| . (2.19) 
The pressure is obtained by 

p = n 2 U -T^ In [l - c -( fc )/ T ] . (2.20) 

k 

For the entropy and energy densities we have the relations 



s 

k 



E 



r ^iVg-lnfl-e- g ( fe )/ T 



(2.21) 



e = Un 2 + t(k)N~ t . (2.22) 



2.5. Accuracy control and fluctuation region 

A necessary condition for the above-outlined relations to apply is the smallncss of the 
parameter 

{V na 3 (d = 3), 
mU (d = 2), (2.23) 
^mlJ/n (d=i). 

At T ^ nil, the parameter 70 directly controls the systematic error of the theory, 
and the condition 70 <C 1 is sufficient for the theory to be accurate. At higher 
temperatures, the condition (|2.23[) is only necessary, since there exists the so-called 
fluctuation region where the fluctuations of the order parameter are essentially non- 
linear and are not captured by our theory. The closeness to the fluctuation region is 
described by the dimensionless parameter 

x = — — P V ( . 2.24) 

{m d T 2 U 2 )—* 

In 2D and 3D systems, fie (T) is the critical value of the chemical potential for a 
given temperature (for the explicit expressions, see section [5]), in ID systems, where 
there is no finite-temperature phase transition (see subsection 15.41 for more detail), 
[li 1 ^ = 0. The theory applies as long as |x| » 1, getting progressively less accurate 
with decreasing \x\. At \x\ < 1, the theory fails to properly describe condensate and 
superfluid densities the values of which being defined by the fluctuating classical-field 
order parameter. The description of other thermodynamic quantities is better since 
the fluctuation contributions to them are of the higher order than the leading (and 
sometimes even sub-leading) terms. 

The estimates for systematic errors for major thermodynamic quantities away 
from and within the fluctuation region are given in sections [5] (in the superfluid region) 
and 17.31 fin the normal region). 
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In the fluctuation region and its vicinity, the theory can be improved by 
incorporating accurate description of fluctuation contributions by — universal to all 
weakly interacting U(l) models and available in literature dimcnsionlcss scaling 
functions. The description is outlined in section [5] 



3. Beliaev diagrammatic technique revisited 

3.1. The Hamiltonian and approach 

The diagrammatic technique for bosons can be derived in terms of functional integrals 
over the classical complex- valued fields ip{v, r). propagating in the imaginary time from 
t = to t = f3 = 1/T (we set h = 1), and subject to the /3-periodicity constraint 
with respect to the variable r [5] . The classical-field grand-canonical Hamiltonian has 
a form 

H = H + H int + H x , (3.1) 

where 

H ° = i / lw ' |2dr ' (3 - 2) 

is the ideal-gas term, with m the particle mass, and 

Hint = \J U{v x - r 2 )|V(r 1 )| 2 |V(r 2 )| 2 dr 1 dr 2 , (3.3) 

is the pairwise interaction term. The H\ term contains linear and bilinear terms 
associated with the chemical potential, /i, external potential, V(r), and the field 7?(r) 
which explicitly breaks the U(l) symmetry, 

iJi = f[V(r) -^]|7/fdr - f [rj* if) + c.c] dr. (3.4) 

For a homogeneous system V(r) = 0. Strictly speaking, at finite rj one cannot refer 
to [i as a chemical potential because the total number of particles is not conserved. 
However, we employ 77-terms solely for explicitly breaking the symmetry and stabilizing 
supercurrent states; at the end of the calculation we take the limit — * 0. We will 
therefore continue referring to [i as the chemical potential. 

Standard finite-temperature treatments of the WIBG typically suffer from infra- 
red divergencies [HI [HI [20l [26j, especially in lower dimensions. The ultra-violet 
divergence is removed by introducing the notion of the pseudopotcntial which allows 
one to express final answers in terms of the scattering length alone. Our treatment 
successfully and systematically deals with divergence issues, and features (at least) 
three important considerations: 

- First, the derivation is based on the appropriate definition of the pseudopotcntial 
U(k) ~ U(0) for low momenta [see the discussion of equation (|3.44[) below] in any 
dimension. 

- Second, the U(l)-symmetry breaking field rj introduces a gap in the spectrum 
of Goldstone modes in a well-controlled way. This suppresses long-range phase 
fluctuations of the order parameter and removes infra-red divergencies in the behaviour 
of correlation functions and self-energies without modifying any of the physically 
important quantities in the relevant order of approximation. In lower dimensions, 
finite rj results in a finite value of the genuine condensate density which dramatically 
simplifies the description. In this respect, a small rj does essentially the same job as 
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V(r)-fi 



Go(r-r') * 



ry*(r) ry(r) 



^in(r) Vout(r) W(r-r') 



Figure 2. Graphical objects representing the single particle propagator G'°'(r — 
r'), the external field V(r) — /i , the symmetry breaking fields ??*(r) and rj(r), 
the condensate lines i/!; n (r) and i/> ou t(r), and the interaction W(r — r') from left 
to right. Here, and in all other figures in the paper, we will only write down the 
coordinates or the momenta of the diagrammatic elements, omitting the frequency 
or imaginary time dependence. 



Berezinskii's finite-size trick, cf. [2T], or Popov's special (hydrodynamic) treatment of 
long- wave parts of the fields [5]. In all final answers one can set r\ = 0. 
- Third, the effects of both interaction and chemical/external potential are 
fundamentally non-pcrturbative separately below the critical temperature. In the 
absence of interactions a positive chemical potential immediately leads to the 
instability of the ideal Bose gas, so that there is no way of consistently introducing a 
positive [i — crucial for describing the system below the critical point — before switching 
on the interaction. This observation explains the idea behind omitting all terms in 
(|3.4[) from the non- interacting Hamiltonian, Hq. At the end of the day, we will use 
fj, — 2nU (0) and T as two variables describing the thermodynamic ensemble. 

3.2. Diagrammatic expansion 

Since the non-interacting Hamiltonian corresponds to the ideal gas at chemical 
potential approaching zero from below, i.e. with no Bose-Einstein condensate, the 
corresponding (bare) Matsubara Green's function reads 



where e(fc) is the single particle dispersion relation. We typically assume a parabolic 
dispersion relation e(k) = k 2 /2m, but our final answers (excluding formulae for 
supcrfluid properties which explicitly invoke Galilean invariance) are valid for arbitrary 
e(k) with a parabolic dependence on momentum in the long-wave limit, e.g. for the 
tight-binding spectrum. We use the Matsubara imaginary-frequency representation: 
£ = £ s = 2sirT (s = 0,±1,±2, ...) with T the temperature. For graphical 
representation of diagrammatic expansions and relations we introduce a set of objects 
in figure depicting the single-particle propagator, the condensate, and various terms 
in the Hamiltonian. 

The non-perturbative response to H luti and H\ is accounted for by Dyson 
summation. First, we consider diagrams that feature only one incoming or outgoing 
line — we call them tails. Given our starting point with no condensate in the non- 
perturbed system such particle-number-changing diagrams exist only due to the 
symmetry breaking field 77. The frequency of the line connecting 77 to the rest of 
the diagram is zero, by frequency conservation. The Dyson summation of all tails 
attached to a given point replaces them with a single line which we denote as V'in(r) 




(3.5) 



G(°)(C,k) = [i£-e(fc)j 



-1 



(3.6) 
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>/w(r) ^ r r > * r r > i c5> 



r r' 



e* ut (r') 



■> — • — — + y^X— + 

Figure 3. Diagrammatic expansion for 113.71 1 and 1 13.81 1. 



r r 



r' r 



( Qin(r') ^- 



and V'out(r), if the tails are incoming and outgoing, respectively, see figure [3] Below 
we write expressions in the frequency representation; if the frequency argument is not 
mentioned it is implied that its value is zero. We also adopt a convention of integration 
over repeated coordinate/momcntum/frequency arguments. The Dyson equation for 
^in(r) then reads: 

^ in (r) = -G (0) (r-r')77(r') + 

G(°)(r - r')[V(v') - /#in(r') + G (0) (r - r')6 in (r') . (3.7) 

Here Q ln is the sum of all other diagrammatic elements attached to the first line which 
are not accounted for by the first two terms, i.e. excluding diagrams with the field n 
and diagrams connected to the first solid line by the [V(r) — fi] vertex. The subscript 
'in' reminds that Oj n has an extra incoming particle line. Similarly, 

Vw(r) = -^(r')G (0) (r'-r) + 

Vw(r')mr') - m]G ( V - r) + 9 out (r')G (t V - r) . (3.8) 

The fact that (r) is a real even function of its argument implies that 

Vw = ip* n , 6 ou t = 0* n ■ (3.9) 

The diagrammatic expansion for the tail is identical to that for the condensate 
wave function defined as an anomalous average 

^o(r) = (V(r)) ee Vin(r) . (3.10) 

Correspondingly, the condensate density is defined as no(r) = | V-'o ( r ) 1 2 - From now on 
we will write ipo for ip- m and -0o f° r V'outj and use the notions of condensate lines and 
tails on equal footing. Of course, in the limit of rj — > speaking of the condensate 
wave-function is meaningful only when the long-range fluctuations of the phase are 
small. 



3.3. Normal and anomalous propagators 

The next step is to introduce exact, or 'bold', particle propagators and work with 
skeleton diagrams. There are a couple of differences between our approach and the 
standard Beliaev technique. Our bare (thin-line) propagators have zero chemical 
potential and this, according to (|3.7|) . immediately results in a constraint relating 
the condensate density to the chemical potential and Oi n [see (|3.21[) below]. Also, the 
chemical/external potential have to be explicitly introduced into the standard Beliaev- 
Dyson equations for the normal and anomalous Green's functions, see e.g., [371 12H1 • 
These equations can be defined purely diagrammatically. To this end — proceeding in 
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Figure 4. Symbols used for normal, G, and anomalous, F, propagators in a 
homogeneous system. 
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Figure 5. Bcliaev-Dyson equations for a homogeneous system. 



the frequency representation for the sake of definitcness — we introduce the normal 
Green's function, G(£,ri,r 2 ), defined as the sum — with a global minus sign, which is 
a mere convention — of all diagrams which have an incoming G^'-line with frequency 
£ to point i"i and an outgoing G^-linc (with the same frequency) from point r 2 . The 
anomalous Green's function, -Fj n (£, ri, r 2 ), by definition, has an incoming 
with frequency £ to point ri and another incoming G^-line with frequency — £, by 
conservation of frequency, to point r 2 . The anomalous Green's function -F ut(£, ri, 1*2) 
is a counterpart of the function -Fi n (£, r 2 ): instead of two incoming it has two 
outgoing -lines; one with frequency £ from point ri and another with frequency 
— £ from point r 2 . The symmetry under exchanging the end points of the anomalous 
Green's functions immediately implies the following relations 



F in (t,r u r 2 ) = F in (-£,r 2 ,ri) , 
-Fout(£,ri,r 2 ) = F out (-£,r 2 ,ri) . 



(3.11) 
(3.12) 



Since complex conjugation is equivalent to changing the sign of the Matsubara 
frequency and direction of propagation we also have 

[G(£,n,r 2 )r = G(-£,r 2 ,n) , (3.13) 

^*(e J ri,r 2 )=F out (-C 1 r 2 ,r 1 )=F out (e,r 1 ,r 2 ) . (3.14) 
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The physical meaning of G, F in , and F out follows from the structure of the two- 
point correlation functions in the imaginary-time-coordinate representation: 

(#ri,ri)^*(75,r a )> =-G(ti-75, n,r 2 ) +^(ri)^(r 2 ), (3.15) 

(^(Ti,r 1 )ip(T2 7 r 2 )) = F in (r! - r 2 , r l5 r 2 ) + Vo(i"i)V , o(i - 2)- (3.16) 

These relations can be readily checked by expanding the averages into diagrammatic 
series. The special case of (|3.15[) . corresponding to T\ — > r 2 — and = r 2 , relates 
the local density to the normal Green's function and the condensate density: 

n(r) = -G(t = -0,r,r) + n<,(r) . (3.17) 

The Beliaev-Dyson equations then read 

G&ri,r a ) = G(°)(e,r 1; r 2 )+ G{^v u v')[V {v') - ^ {^v\v 2 ) 

+G(C,r 1 ,r')S 11 (e,r',r")GW(e,r",r 2 ) 

+ F in (£, rr, r') £ 20 (£, r', r") G^fo r", r 2 ) , (3.18) 

Fin^ri.ra) = F in (£, n, r')[^(r') - M ]G (0) (-£, r 2 , r') 
+ F in (£, r l5 r') E u (& r', r") GC >(-£, r 2 , r") 

+ G(£, r l5 r') E 02 (£, r', r") G (0) (-£, r 2 , r") , (3.19) 

with the standard definition of self-energies E's as sums of diagrams which can not be 
cut through a single G or F line. Equations (|3.18[) and (|3.19|) are shown graphically 
in figure [5] for a homogeneous system in momentum representation. The complexity 
of the theoretical solution is in the evaluation of the and E functions. 

3. 4- The chemical potential and the Hug enholtz- Pines relation 

Since the bare Green's function at zero frequency is identical to the inverse Laplacian 
operator one can cast (|3.7|) in the differential form 

- M*) + [V(r) - fi] Mr) + e in (r) = ry(r) . (3.20) 
2m 

This equation reduces to the Gross-Pitaevskii equation at low enough temperature 
when the leading term in 0i n is cx |V'o( r )| 2 V'o( r )- In the homogeneous case at rj — > 0, 
when ipo and 0j n are coordinate- independent, equation (|3.20p simplifies to 

n = e in /v^o = e O ut/^ *: (3-2i) 

generalizing relations (21.3)-(21.5), and fi = (0\d(H in t/V)/dn o \0) , in [2SJ to finite 
temperatures. 

There is also an exact relation between En, Eo 2 , 6i n , and ipo (see also [3]). 
Here (and only here!) we assume that all diagrams for E's are in terms of G^ and 
condensate lines. Let be the sum of diagrams contributing to 0j n with I incoming 
and / — 1 outgoing condensate lines. Then (for £ = 0) 

oo 

E n (r,r')^o(r') = ]T ID® , (3.22) 

2 = 1 

because each diagram, with I — 1 incoming condensate lines, contributing to Sn(r, r') 
produces — upon integration over r' with the weight V'o( r ') — a diagram contributing 
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to , and there are I such diagrams contributing to E u . An identical argument 
leads to 

oc 

£o 2 (r,r')VoV) =J2V-V D $ ■ ( 3 - 23 ) 

1=2 

By subtracting (|3.23[) from (|3.22[) we obtain 

En(r,r') Vo(r') - £ 02 (r, r') VoV) = e in (r) . (3.24) 

In the homogeneous case ^>o = V'o- Then, in the 77 — * limit, equations (|3. 241) 
and (|3.2ip can be combined to yield the Hugcnholtz-Pines relation [4] (and its finite- 
temperature version [5]) 

H = £ u (0,0)- £02(0,0) . (3.25) 



3.5. Beliaev-Dyson equations in the presence of homogeneous superflow 

In order to discuss the superfluid properties of the homogeneous system, we add a 
phase factor to the symmetry breaking field 

r ? (r)=7 ? oe iko ' r . (3.26) 

which readily translates into phase of the condensate wave function 

Mr) = V^e ik - r , (3.27) 

The only difference with the previous discussion is that now we have to associate 
a finite momentum ±ko carried by the condensate lines and modify the momentum 
conservation laws accordingly. Then (|3 . 2 1 [) and (|3.25|) become 

fx = 6/Vno" + fco/2m - WV"o > (3-28) 
£ 11 (k )-£o2(0) = ^= = ^-^ + ^, (3-29) 



where 9 in = 6 ou t = ©• 

It is convenient (for transparency of expressions which follow) to combine 
frequency and momentum into a single "4-momentum" variable P = (^, k) and to 
introduce an auxiliary momentum P' = (— ^,2ko — k). The symmetry between the 
two ends of the anomalous Green's functions and, equivalently, between the two ends 
of the anomalous self-energies is then expressed by (accounting for the momentum 
carried by the condensate lines) 

Fin/out (P) = Fin/out (P') , (3-30) 

£20/02 (P) - £20/02 (-P') ■ (3.31) 

[In a more comprehensive notation scheme one has to mention momenta of both 
incoming lines in F in/ont (P, P').] 

Complex conjugation of propagators and condensate lines changes the signs of 
their 4-momcnta. This property can be used to prove the symmetry relation for the 
Green's function 

G*(P)| k0 =G(-P)U . (3.32) 

Similar symmetry relations take place for the anomalous Green's functions and all 
three self-energies. 
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In the momentum representation, inverting the direction of k does not change 
the analytical expression for propagators. Hence, by inverting the direction of all the 
lines, including the condensate ones, it follows that 

Fin(P) = F out (P) = F(P) , (3.33) 

£ 20 (P) - £ 02 (P) ee £(P) . (3.34) 

(For the normal Green's function, G(P), and the self-energy, £ 1:L (P), inverting 
momentum directions results in the same series, and in this sense is trivial.) 

We are ready to formulate the pair of Beliaev-Dyson equations in the momentum 
representation. Using symmetry properties into consideration and the shorthand 
notation S = Su we obtain 

G{P) = G (0) (P) + G(P)[E(P) - ^]G (0) (P) 

+F(P)£(P)G (0) (P), (3.35) 
F(P) = P(P)[S(P')-/i]G (0) (P') + G(P)£(P)G {0) (P') . (3.36) 
The solution in terms of self-energies reads 

G(P) = !i±fajj^WHjf , ,3.37) 

FlF) = -§m' ( " s) 

where 

D = E 2 (P)-[e(\2k -k\) + ^(P')- f i+^}[e(k) + E(P)-fi-^}. (3.39) 

With these relations at hand, one can calculate the current density induced by the 
phase gradient in the condensate wave- function, see subsection 14.51 



3.6. Low-density limit in 3D and 2D: Pseudo-potentials and scattering lengths 

In two and three dimensions, the expansion in terms of the bare interaction potential 
can be (and in most realistic cases, is) non-perturbative. The system is regarded as 
weakly-interacting only because of the low density of particles. In one dimension, 
the physics is perturbativc instead in the high-density limit for a given potential. 
Theoretically, dealing with the strong bare potential implies summation of an infinite 
sequence of ladder diagrams and produces an effective interaction in the form of the 
four-point vertex [2], I\ see figure [6j The analytical relation behind figure [6] reads 

r(p 1 ,p 2 ,Q) = 

U(q) + w (q - k ) G ( p i + K ) G{P 2 -K) T(P U P 2 ,K) (3.40) 

K 

ee U(q)+J2 T(P 1 +K,P 2 -K,Q - K)G(Px+K)G(P 2 -K)U(k) . 

K 

When the bare-interaction lines are replaced with F's, the rest of the series becomes 
perturbativc (excluding the critical region). On the technical side, working with 
T's is convenient as long as one does not intend to systematically take into account 
higher-order corrections, utilizing thus only the (simple and transparent) leading-order 
expression for T(0,0, 0). For the higher-order corrections, equation (|3.40[) involves 
three different length scales: the size of the potential, Rq, the healing length of the Bose 
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Figure 6. Ladder diagrams leading to equation 113. 401 1. 



condensed system, and the de Broglie wavelength, while only the non-perturbative 
physics at the scale Rq requires summation of the ladder diagrams. Hence, it makes 
sense to construct an object with a much simpler structure than T that captures the 
non-perturbative physics (and thus coincides with the leading-order expression for 
r(0, 0, 0)), and then systematically investigate the difference between this object and 
r. To achieve this goal let us define the pseudo-pontential J7(q) by the equation (see 
figure [7]) 

U(q) = W(q) + J2 W ^ - k ) n ( k ) U Q4 ( 3 - 41 ) 
k 

= W(q)+£ U"(q-k)n(k)M(k) , 

k 

where IT(k) is such that when k is much larger than the inverse healing length or 
thermal momentum the value of n(k) approaches that of the product G(P\+K)G(P2— 
K) summed over the frequency of the 4-vcctor K. That is, 

n(k) -► -t4tt at k -> oo. (3.42) 

In this case, T(Pi, Q) ~ U (q), and one can expand the difference in a perturbative 
series. As a result, we arrive at the diagrammatic technique where instead of thin 
dashed lines standing for the bare interaction potential we have bold dashed lines 
representing the pseudo-potential U, with an additional requirement that whenever 
two (normal) Green's function lines are sandwiched between two pseudo-potential 
lines, the former are supposed to be summed over the frequency difference (which is 
always possible in view of the frequency- independence of the pseudo-potential), and 
then II has to be subtracted from the result of summation. Specifically, if Pi and P2 
are the two external 4- momenta of the above mentioned diagrammatic element, then 
the following replacement is supposed to take place for internal propagator lines 

Y J G{Pi+K)G(P 2 -K) ^G(Pi+if)G(P 2 -X)-n(k), (3.43) 

where £^ K ' is the frequency of the 4-momentum K. 

As long as we are not interested in the non-universal ultraviolet corrections, we 
can replace ?7(q) with U(0), the systematic error introduced by the replacement being 
controlled by the following dimensional estimate 

C/(q) = U(0) [1 + 0{q 2 Rl)] . (3.44) 

One may wonder how this estimate is reconciled with the momentum dependence of 
T, which is of the order qR^ in 3D, and of the order 1/ \n{l/qR ) in 2D. The solution 
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Figure 7. The diagrammatic expression for the pseudo-potential (left hand side) 
involves the bare potential (first term on the right hand side) and a modified 
two-particle propagator II (k). 



is provided by (|3.43[) explaining that this dependence, which is both universal and 
perturbative, is taken into account by the second-order ladder-type diagram in U. 
In 3D, a natural choice for II (k) is 

m = ~ (* = 3), (3.45) 

in which case we have the usual [23 [55] (see also below) 

, , 47ra , . 

U(0) = (d = 3) , (3.46) 

m 

with a the s-wave scattering length. In 2D one cannot use (|3.45|) because of the 
infra-red divergence of the integral in (|3.41[) . A reasonable choice here is 

n(k) = -\ 1 (d = 2), (3.47) 

2 e(k) + e 

where eo is an arbitrary low-energy cutoff. The particular value of eo is not important 
since final answers are not sensitive to it. Within the systematic-error bars of 
the pseudo-Hamiltonian description discussed below, any choice of e such that 
nil < eo ^ n/m is equally reasonable in terms of accuracy, provided the temperature 
is not much larger than n/m; moreover, replacing eo with ck is also acceptable. (An 
optimal choice of e in the regime T ^> n/m will be discussed in subsection l7.2n There 
is also no need to introduce II in d — 1. We will assume that formally II(k) = in 
d = 1 in which case our final answers can be used as written in all spatial dimensions. 

Given that (|3.45[) and ()3.47j) are not unique [there is a free parameter in l|3.47[) ]. 
it is instructive to explicitly relate two pseudo-potentials, U± and U2 — corresponding 
to ITi and 1T2, respectively — to each other. We notice that (|3.41|) implies 

fa(q) = f/i(q)+^f/i(q-k)[n 2 (k)-n 1 (k)]C/ 2 (k). (3.48) 

k 

In view of (|3.44p this simplifies [up to 0(q 2 B%) terms that we neglect in what follows] 
to 

U 2 = Ui + U X C X2 U 2 , (3.49) 
C 12 =^[n 2 (k)-n 1 (k)]. (3.50) 



k 



If 111 and II2 are defined by (|3.47|) with different cut-offs e ' and e , we have 

(2) 

TYl f 

Ci2 = ^ln^y (d = 2). (3.51) 
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The right choice of the functions IT and n 2 implies that U\ and U2 differ only by 
sub-leading terms. Thus, we can expand the r.h.s. of (|3.49[) in powers of | C/i C12 1 <C 1: 

U 2 = Ui + UlC l2 + ... . (3.52) 

Formally, the expansion in terms of the pseudo-potential is perturbative only as 
long as we exclude the high-momentum contributions to the (M > 3)-body diagrams 
generating, upon complete summation, M-body scattering amplitudes. In a dilute 
gas, the corresponding diagrams are small (contain extra powers of the gas parameter 
na ) and are neglected in this manuscript. 

The expression (|3.46[) has the meaning of mapping a dilute three-dimensional 
system with an arbitrary short-range interaction potential onto a system with the 
interaction potential (|2.1|) . so that the pseudo-potentials of the two systems coincide. 
The same approach is possible (and popular) in 2D, the parameter a being called the 
two-dimensional scattering length, since the mapping applies to scattering properties 
as well. For a given potential hi, the value of a can be obtained either from the 
asymptotic behaviour of the pseudo-potential U at appropriately small eo, or from the 
asymptotic behaviour of the scattering amplitude /(k, k') at fc, fc' — > [28]. Both 
asymptotic behaviours are closely related to each other since the large-fc limit of the 
kernel II(k) in (|3.41[) coincides with the large-fc limit of the scattering kernel 

m 

lW)= fc , 2 _ fc2 + iQ (3-53) 

in the integral equation for /(k, k'), 

/(k, k') = W(k-k') + W(k-q)n sc ( g , fc') /(q, k') (3.54) 
q 

= W(k-k') + ]T f(k, q)n sc (q, fc') W(q-k') . 
q 

Moreover, a direct relationship between U and /(k, k') is readily obtained by noticing 
that lESjUl and (|5^T|) imply [cf. (338J)] 

/(k,k') = C/(k-k')+^(7(k-q)[n sc (< 7 ,fc')-n( ( z)]/(q,k'), (3.55) 
q 

which dramatically simplifies upon replacing t/(q) with U(0), in accordance with 
(|3.44p . In this case, /(k, k') is k- independent, and for / = f(k') we get 

/(fc') = u + uf(k') [n sc (<7, fc') - n( q )] . (3.56) 

q 

Substituting (|3.47[) and p.53p into (|3.56p and performing the integral, we find 

f(k ') = V?Z? f3 57) 

ln( v / 2me r ;/fc') +27r/mU' + i7r/2 ' v ' ' 

Comparing this to the known hard-disk result (see, e.g., [29] ) 

/(fc / ) = 1 /9 , (3-58) 

where 7 = 0.5772 ... is Euler's constant, we conclude that 

U= - , o/ f /m - . (3.59) 
ln(2/a 2 me ) -27 y J 
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Figure 8. The lowest order diagrams for the self-energy S in 114. H and the 
anomolous seflenergy S in 114.21 1. 



In 3D, equation (|3.56|) — with II(k) given by (|3.45|) — yields 

thus leading to (|3.46|) by comparison with known expression for the hard-sphere 
scattering amplitude. 

4. Thermodynamic functions in the low-temperature region 

4-1- Basic relations and notions 

Explicitly calculating the lowest-order diagrams shown in figure [8] and utilizing (|3. 17j) 
one finds 

£(P) = -2G(r = 0, t = -Q)U + 2n Q U = 2nU , (4.1) 

E(P) = n U . (4.2) 

We see that within the first approximation, both £ and E turn out to be momentum- 
and frequency-independent. [It is easy to check that the next-order diagrams 
inevitably introduce momentum and frequency dependence to self-energies and 
drastically change the structure of the theory.] At this level of accuracy, the chemical 
potential equals to fx = 2nU — n U according to the Hugenholtz-Pines relation. As 
mentioned above, it is extremely convenient to use an effective chemical potential 

(1 = (i- 2nU , (4.3) 

as a thermodynamic variable to characterize properties of the WIBG. [Note that jl is 
negative.] Within the same accuracy we can substitute E with —jl = \jl\ and simplify 
expressions for G and F to 

F ^ = wrlm- (4 - 5) 

£ 2 (k) = e(t)[ e (t)+2|A|]. (46) 
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A note is in order here. Below we will calculate higher-order corrections to the 
chemical potential which are necessary for the construction of the accurate equation of 
state, see (|4.9p and (|4.11l) . However, it is still possible to use fx in the definitions of the 
propagators and the spectrum of elementary excitations while keeping the accuracy of 
the entire scheme intact, see, subsection [5] where we study systematic errors involved 
in approximations. Here we mention briefly that the anomalous average contribution 
to thermodynamic properties is always small, and thus further corrections to F are 
negligible. The same is true for G at low temperature. At temperatures T 3> noli, 
on the other hand, one can ignore tiny modifications in the spectrum because of an 
additional small parameter uqU/T. 

For the total density we have 

n = n - G{r = 0,r = -0) 

= * + E[1|ifV + ^)-i], (4.7) 

where 



N e = U' T - l) (4.8) 



is the Bose-Einstein occupation number for the mode with the energy e. Formula (|4.7p 
includes the leading and sub-leading terms. To get an expression for the chemical 
potential with the same degree of accuracy, we need to take into account higher-order 
diagrams and go beyond the leading-order expressions for S and E. This is because 
in the absence of interaction the chemical potential is identically zero and thus the 
leading term is proportional to U. 

There are three second-order diagrams contributing to 0. The first one is 
the anomalous Green's function convoluted with the bare interaction potential (see 
figure O The second one is the 'sunrise' diagram with the proper correction (|3.43D 
for the ladder structure and the third one is similar to the sunrise diagram but with 
propagators connecting different vertices (see the two diagrams in the second row in 
figure \§§ . The latter two can be safely neglected away from the fluctuation regime 
because they involve an additional small parameter 70 or 7y, see Section [S] for the 
analysis and definitions of the fluctuation region and parameters 7 in different spatial 
dimensions. It is worth noting that consistently taking into account contributions 
of the two neglected diagrams in the condensed regime would require simultaneously 
going to higher orders in the expansions for self-energies, figured] The latter, however, 
is impossible without sacrificing the attractive paradigm of independent quasiparticles 
with Bogoliubov dispersion. 

Keeping the leading and the largest sub- leading diagrams for 0, we get 

/i = -2G(T = -0,r = 0)[/ + n W(0)-^W(q)F(q,T = 0). (4.9) 

q 

We have no choice but to use the bare potential here because all ladder diagrams 
leading to the pseudopotential vertex are already absorbed in the anomalous Green's 
function, by construction of the latter. This feels unsatisfactory only at first glance 
since simple formal manipulations allow one to express (|4.9|) in terms of U alone. 
Let us introduce an auxiliary function A(q) defined by the integral equation (shown 
graphically in figure [TO]) 

A(q) = F(q,r = 0)+n II(q)C/(q)-^ II(q-k)tf(q-k)A(k) , (4.10) 
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Figure 9. Diagrams contributing to the chemical potential fj, of a condensed gas, 
up to the second order. The last two diagrams [that actually have to be corrected 
according to H3.43H are smaller than the previous ones and will be neglected; see 
the text for more detail. 



Figure 10. Definition of A(q). 



It can be used, in combination with the definition of the pseudopotential Eq. (|3.41D , 
to transform the last two terms in (|4.9[) 

n o W(0)-^W(q)F(q,T = 0) = 
q 

n o U(0) - Y, ^(q)A(q) w (no - A Q )U , (4.11) 
q 

where 

A = Y A(q) . (4.12) 
q 

The last approximate equality in (|4.11[) comes from U(q) w U(0) and the observation 
that A(q) vanishes at momenta much smaller than I/Rq. 

The smallness of the parameter |[/IIfc d | at momenta k -C l/-Ro allows one to 
expand A [see (|4.10[) ] in the series: 

A(q) = [F(q, r = 0) + n n(q)C7(q)] + . . . . (4.13) 
For the effective chemical potential we thus get 

fi = -(n + A )U , (4.14) 
where within our order of accuracy we can take 

A o = T = °) + n oCHI(q)] . (4.15) 

q 

We find it convenient to introduce a quantity with the dimension of density, 

n* = -Jx/U = n + A . (4.16) 
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With this quantity the form of certain thermodynamic relations simplifies. For 
example, using (|4.16[) we can replace Uq with n* in (|4.7[) and arrive at the result 

n = n*+n' , (4.17) 



n = 

k 



E 



e{k) - E{k) > + 



2£(fc) ^ v ' E(k) 



(4.18) 



This completes the self-consistent theory because we obtain a closed set of relations 
which define n = n(/x, T) in the parametric form — given some /2, or n*, one calculates 
n from (|4.17|) and (|4. 18)) and then determines fi from (|4.3[) . The integral in (|4.18p is 
convergent not only in 3D, but also in 2D and ID. Hence, at this point the quantity 
rj can be set equal to zero. 

We want to emphasize the fact that all specific expressions derived in this section 
feature a two-parametric representational freedom, within the same order of accuracy. 
First, it is possible to use p with or without higher-order corrections in the spectrum 
of elementary excitations and propagators. Second, there is a freedom of choosing the 
function II(fc). For example, one could be tempted to absorb the first two terms in 
the integral in (|4.18p into the definition of U, such that n'(T = 0) gets identically 
equal to zero, and n^T = 0) equals to the total number of particles. We, however, do 
not see any merit in this protocol, because U becomes dependent on p and cannot be 
considered as a fixed external parameter. Finally, this and analogous "improvements" 
do not make the theory more/less accurate since the difference is of the same order as 
omitted diagrams. 

A remark is in order here concerning the two-dimensional case, where the value 
of U cannot be defined irrespectively of the system density. Even in this case, we can 
proceed with formally independent parameter eo in (|3.47[) till we arrive at the final 
answers along with the estimates of neglected higher-order terms (the latter being 
essentially eo-dependent). After that, the value of eo is selected in such a way that 
the order-of-magnitudc values of neglected terms are minimal. 



4-2. Pressure 



To derive relations for other thermodynamic properties, we start with the pressure as 
a function of T and p. Using the general thermodynamic formula 

_ dp(ji,T) 



Of i 



(4.19) 



and adopting — throughout this subsection — the convention that T is treated as a fixed 
constant, so that partial derivatives with respect to either /j, or p can be replaced with 
ordinary ones, we write 



P = Pc + 



" d» 
n — d/i, 



where 



nzU - T 



EM 



l-e 



-e(k)/T 



is the value of pressure at the mean-field critical density 



,{k)/T 



(4.20) 
(4.21) 

(4.22) 
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Equations (|4.3[) and (|4. 17)) allow one to represent (|4.20j) as 

p = p c - n 2 U + p 2 /2U - 2n'p + {n'fU + n' dp 



23 



(4.23) 



(after doing straightforward integrations by parts). The integral in (|4.23[) is readily 
done by noticing that 

e(fc) dE(fc) 



E(k) 
e{k)N E 



dp 



d 



and thus 



= -T — In 
E(k) dp 



^ k 



1 _ Q -E{k)/T 



-T 



1 _ e ~E{k)/T 



dp 

where e(k) makes the hrst term convergent. The result of the integration is 



n dp 

E 



{ E(k)-e(k)+p-p*Tl(k) 1 - c- E ^> T 

2 1 1 — e-W/T 



I 

and we finally get 

p = fi 2 /2U - 2n'p+{n') 2 U 

-- { E ( k ) - < k ) + A - P 2 ^{k) + 2Tln 



(4.24) 
(4.25) 



(4.26) 



1-e 



-E(k)/T 



(4.27) 



} • (4-28) 



4-3. Entropy and energy 

Whenever n, [i, and p are specified as functions of (T, x), where a; is a quantity of 
arbitrary nature, the expressions for entropy per unit volume, s, and energy per unit 
volume, £, are readily found from the following two generic thermodynamic relations: 



dp 
df 

e = [in + Ts — p 
With x = p we thus find 

- = E 



dp, 
df 



P 2 



= j^-n'p + n u U 



~ £ [E(k)(2N E + 1) - e(k) +P- p 2 U(k)] 



(4.29) 
(4.30) 

(4.31) 



(4.32) 
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4-4- Explicit integrations. T = case 

It is useful to note that the following two integrals (we use ft = —\£l\ 
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(4.33) 



E{k) 



- 2|/i|n(fc) - i 



si} 



a|A| 



(4.34) 



entering relations for thermodynamic quantities — see subsections 14.1114.31 — can be 
explicitly performed. [Here we assume that in 3D and 2D the kernel H(k) is fixed by 
expressions Q3.45P and (|3.47[) . respectively, and remind that in ID it is zero.] Explicitly 
doing the integrals is especially useful at T = 0. In this case, all other integrals in the 
expressions for thermodynamic functions nullify, and the answers reduce to algebraic 
relations. In the grand-canonical form, these relations read 



(4.35) 
(4.36) 

(4.37) 
M 

in the arguments of I^ 1 and 1% , since the integrals are responsible for sub-leading 
corrections, while the sub-sub-lcading terms are beyond our control. 
Performing the integrations, one finds: 



(T = 0) , 
(T = 0) , 

(T = 0) . 



Here we take into account that within our level of accuracy we can replace \jl\ 



l[ d=1 \ 



r(<*=2), 



r(d=3) 



and differentiating with 



r(d=i) 



Jd=2) 



,(d=3) 



W = 



2v^|/i| 3/2 



3?r 



Ml) 



/'I) 



mjl 2 

~8^r 



In- 



2eo 



8m 3 / 2 |A| 5/2 



15tt 2 

respect to |/2| we find 



f<\) 

m 



47T 



In - 



2e 



4m 3 / 2 |/i| 3/2 



3tt 2 



(4.38) 
(4.39) 

(4.40) 

(4.41) 
(4.42) 

(4.43) 



In 2D, we have a freedom of fine-tuning e to simplify the form of the answer. A 
natural choice, especially convenient for the T = limit, is to set 

e = (e/2)|A|, (4.44) 
in which case we have 



r(d=2) _ 



0, 



{d=2) 



mjl 2 
' 16tt 



(4.45) 
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At T — this translates into very compact grand-canonical expressions :29j 

\ii\=fi, n = £ (d=2, T = 0), (4.46) 

P=^r( 1 + ^) ( d = 2 > r = °)> ( 4 - 4 7) 

with 



2U V 8tt 



^= ^(4/^) -27-1 - (d = 2 ' T " Tc) - (4 ' 48) 
It is seen that the simplicity of the form of (|4.46[) comes at the expense of more 
sophisticated (fine-tuned) form of the effective interaction (|4.48[) . If, instead, one 
would opt to simplify the form for the effective interaction, getting rid of sub- 
logarithmic terms: U = (47r/m)/ln(l/a 2 m/x), then (|4.46[) would acquire the generic 
form (|4.35p - (|4.36p . Needless to note that this does not change the sum of leading and 
sub-leading terms in the equations of state since, by construction, the result cannot 
depend on the specific choice of eo- 

4-5. Superfluid density 

The standard way of calculating the superfluid density within the quasi-particlc picture 
is based on Landau's formula for the normal component density, n n , and the relation 
n s = n — n n . Here we employ a more general approach based on the current induced 
by the gradient of the condensate wave-function. By definition, the superfluid density 
is the linear-response coefficient relating the persistent-current density to the gradient 
of the phase, ip, of the (complex) order parameter field: 

j = -V^ = n>. (4.49) 

TO TO 

In the last equality we assumed that <p = ko • r. On the other hand, the average 
current can be calculated from the microscopic operator expression in terms of the 
condensate density and the Green's function 

j = - J- ($t V ^ - H.c.]\ = ^n Q + V - jV(k , k), (4.50) 

2TO \ /TO TO 

k 

where 

7V(k,k ) = G(t = -0,k)| ko . (4.51) 
In the limit of ko — > we obtain the required relation 

n s = n Q + lim V 7V(k , k) . (4.52) 

k u 

At first glance equation (|4.52j) docs not resemble Landau's formula at all, since 
the first term is given by the condensate density, not n. However, the structure of the 
normal average contribution is such that it has a part which completes no to the total 
density and a part which equals (with minus sign) to the normal density component. 
To derive this result, we start with the expression for the Green's function 

n( t i£ + e(k-2k )-/2 

k k = — — -— , 4.53) 

M + K + e(k-2k ) - nM - 6 k ) + fj] 



Beliaev technique for a weakly interacting Bose gas 



26 



and solve for the roots of the denominator in (|4.53[) 

e(fc)- e (k-2k ) ^ 
Ki, 2 = ± 



'e(k)e(k-2ko) - [e(k)+e(k-2ko)]/i ■ 



e(k) - e(k-2k ) 



to rewrite the Green's function identically as 
a 1 — a 

Gr 



i£ -Ri i£-R2 
e(k - 2k ) - /i - Ri 



(4.54) 

(4.55) 
(4.56) 



Rl — i?2 

Now we express frequency sums through the Bose functions to arrive at 

iV(k)| ko = aAT(-Ri) + (a - 1)(1 + N(-R 2 )) . (4.57) 

The rest of the calculation is straightforward. In view of the limit to be taken in 
(|4.52|) it is sufficient to keep only terms linear in [e(k) — e(k— 2ko)] ~ 2k • ko/m in 
Ri,2, ct, and iV(k)|k . We omit the algebra of doing the expansion which leads to the 
result 

ON 



n 



' J rim 



/< 2 



dm E 2 (k) 



1 + 2N E _ 
2E{k) dE{k) 



(4.58) 



Finally, using (|4. 7[) to exclude uq and integrating by parts to simplify the expression 
we arrive at Landau's formula 

— ■ (4.59) 



ly, k 2 dN 
m~dE ' 

k 



5. Expansion parameters. Estimates for higher-order terms 

Let us analyze the structure of small parameters which control the applicability 
of the diagrammatic technique presented above. As we will see a posteriori, it is 
sufficient to consider two characteristic limits: (i) the T ~ case and (ii) the finite- 
T contributions of diagrams with zero frequencies (i.e., classical-field contributions). 
The analysis is based on dimensional estimates of the diagrams, the fact that an extra 
interaction vertex increases the total number of propagators by two, and that the 
largest contributions come from small momenta (and frequencies), where the normal 
and anomalous propagators behave as 

\G\ * \F\ « e -j^L , c 2 = M/m, (5.1) 

and we do not need to distinguish between them. We will further assume that there 
is an infra-red momentum cutoff k\ <C kg where 

fc = VrrmU (5.2) 

is directly related to the inverse healing length. The symmetry-breaking field r/ also 
changes the behaviour of propagators at k <§C kg , but for our purposes we do not need 
the explicit form of propagators at finite rj because we will keep k\ large enough to 
neglect effects originating from finite rj. 
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At T < T c there are two kinds of generic vertices in higher-order diagrams, 
namely (i) full vertices, where four propagators meet, and (ii) vertices with one 
condensate line (vertices with two and three condensate lines are all accounted for by 
the lowest-order diagrams for the self-energies and can not be part of the diagrammatic 
expansion) . A naive definition of the diagram order as the total number of vertices 
proves inconvenient, since contributions of condensate vertices turn out to be larger 
than contributions of full vertices. Below we will see that the difference is exactly 
compensated by the fact that condensate vertices (generically) come in pairs. This 
suggests to define the diagram order as the sum of the total number of full vertices 
and half the total number of condensate vertices (plus 1/2 for anomalous diagrams). 



5.1. T = case 

At T = 0, the structure for the dimensionlcss factor associated with adding an extra 
full vertex to a diagram is 

7 f U) ~ uj 6(Ak)S(A0 [\G\ d d kdtf , (5.3) 

with S(Ak) and <$(A£) representing the (5-functions taking care of momentum and 
frequency conservation laws. The dimensional estimate, following from (|5.3|) and (|5.1|) . 

reads 

(full) nV*(mU?l* 
fc?- 

The structure for the dimensionlcss factor associated with adding an extra pair 
of condensate vertices is 



i { r ] - rz' ■ (5-4) 



7 < cn } ~ U 2 n J [S(Ak)S(A0] 2 [\G\ d d k d£] , (5.5) 

and from the corresponding dimensional estimate one can see that 

7 ( cnd ) „ 7 ( M ) rnnU _ ^ 



Equations (|5.4[) and (|5.6[) clearly show the infrared problem of the theory in dimensions 
d < 3 [THlUni HO] (as usual, a zero power of the momentum cutoff should be understood 
as a logarithm), and emphasize the usefulness of the cutoff-enforcing field ij. For the 
diagrammatic expansion to be consistent, we need k\ to be large enough in order to 
guarantee 70 < 1. On the other hand, we do not want the field rj to significantly 
affect the physics of the system, and thus need k\ <C ko. Hence, the smallest possible 
expansion parameter one can afford without distorting the physics corresponds to 
ki < ko, in which case 

!V na 3 (d = 3), 

mU (d = 2), (5.7) 
^/mU/n (d=l). 

Here we took into account that \jl(T = 0)| ~ nil, and expressed U in terms of a 
in 3D. It is clear that 70, equation (|5.7j) . is the actual expansion parameter for all 
local thermodynamic quantities. These quantities, by referring to large length scales 
at T = 0, should mainly have contributions from wavevectors < ko. Technically, it 
means that infrared divergencies of individual diagrams have to cancel each other in 
the final answers and the resulting integrals become convergent at the wavevectors 
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k ~ fc , leading to a well defined expansion in powers of 70 (|5.7p . [Since dimensional 
estimates do not distinguish between pure powers and powers with logarithmic pre- 
factors, we extend the meaning of the term 'in powers' to the powers with logarithmic 
pre- factors.] 

Note that all our relations between the thermodynamic quantities that do not 
vanish in the T = limit (specifically, n, fi, p, e, and no), are accurate up to the first- 
order corrections in 70, so that their systematic errors are of the order of 7q (with 
possible logarithmic pre- factors). The relation for the entropy is accurate only up to 
the leading term, which actually has the same quasi-particle origin as the sub-leading 
terms in all non- vanishing at T — > quantities. The same is true for the heat capacity 
that behaves similarly to entropy at T — > 0. 



5.2. Finite-T zero-frequency contributions 

The dimensionless factor associated with adding an extra full vertex to a diagram 
consisting of zero-frequency propagators scales as 

7 £ uU ) ~ (U/T) J S(Ak) [\G\Td d k] 2 , (5.8) 

yielding the dimensional estimate 

For the pair of condensate vertices we have 

7 (cnd) _ {U/T) 2 nQ f [S{Ak)] 2 [IGjrd^] 3 . (5 . 10) 



Similarly to the T = case (assuming that \fl\ ~ n U, we get 

7 (cnd) ^ ^W,) m^U_ (5 n) 

In complete analogy with the T = case, the largest possible k\ (enforced by 77) 
cannot exceed \/mnoU. Assuming that fei < s/mrioU results in 

m / „ \ 2—d/2 
(full) (end) 1 n\ 

nU \n J 

This expression shows that at T <C nU the parameter 70 dominates over 77 1 , and thus 
this temperature regime is equivalent to T = 0. In the crossover regime, T ~ nU, both 
70 and jt are of the same order. In the regime T 3> nU , the classical-field parameter 
7t dominates, and becomes of order unity when no gets sufficiently small. The latter 
situation corresponds to the fluctuation region, where the perturbative theory becomes 
inadequate, since it fails to properly describe the non-linear long-wave fluctuations of 
the classical component of the quantum field. Before hitting the region ■jt ~ 1, one has 
to switch to the description in terms of scaling functions [151 117] which are universal 
for all U(l) weakly-interacting theories, see section [TU1 

At T > nU, all our relations for the thermodynamic quantities are accurate up 
to the first-order corrections in 77 1 , so that their systematic errors are of the order of 
7|, (with possible logarithmic pre- factors). Note that at T ^ nU the leading terms 
for all the thermodynamic quantities, except for the chemical potential (and with a 
reservation for the density in ID where the corresponding condition is T > nU/jo), 
are the same as for the ideal gas. 
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To summarize our analysis, we present below an order-of-magnitudc estimate of 
systematic uncertainties due to omitted higher-order terms for each thermodynamic 
quantity. It is convenient to start with fi and recall that the omitted term in (|4.9[) 
comes from the 'sunrise' diagram for O. A dimensional analysis of this diagram gives 

A/i/m ~ max [7o>7T] ■ (5.13) 
which transforms into a similar estimate for the total and superfluid density 

An/n* ~ max[7o,77-] . (5-14) 

An s ~ An . (5.15) 

This can be seen from the dimensional analysis of the second-order diagrams 
contributing to the self-energies. Relation (|4.20p for the pressure implies 

Ap ~ n/imax[7Q,7^] . (5.16) 

In the case of entropy, it is convenient to start with the regime T > nil when — using 
(|4.29[) to relate the uncertainty in s to Ap, An, and Au — we find 

As ~ ^n 7 ^ (T>nU). (5.17) 

We see that at T ~ nU, in contrast to the behaviour of other thermodynamic functions, 
the uncertainty in s scales only as the first power of 70 (simply because s itself gets 
of order ~ 70). This scaling persists down to T — > 0: 

As/s - 7o (T< nU) , (5.18) 

because, at T <C nU, the thermodynamics of the system corresponds to the 
generic low-temperature behaviour of superfluids, where the leading temperature- 
dependent contributions are due to phonons. Hence, the first-order correction to 
the sound velocity, which is of the order of 70 as is seen from the expression for the 
energy, translates into the order-70 correction to the entropy. At T <C nU, when 
thermodynamics is exhausted by dilute non-interacting phonons, the correction to 
the sound velocity can be found directly from the zero-temperature compressibility, 
and the accuracy of the expression for s (and other thermodynamic quantities) can be 
improved. At T ~ nU, however, the order-70 correction to the quasiparticle dispersion 
law goes beyond the Bogoliubov ansatz [30] . Note also that improving the value of 
the sound velocity in the T <C nU, while rendering the answers for thermodynamic 
quantities more accurate, is inconsistent with retaining the Bogoliubov form of the 
spectrum for all the momenta. 

Finally, equation (|4.32p in combination with the above results yields the estimate 
for the uncertainty of energy: 

Ae ~ Ap ~ n/2max[7o,7f.] . (5.19) 



5.3. Fluctuational contributions 

In all spatial dimensions, there is an interval (in terms of temperature, if density is kept 
fixed, or in terms of chemical potential at fixed temperature, etc.) where jt becomes 
of order unity, and the systematic perturbative description breaks down. In 3D and 
2D this happens in the vicinity of the superfluid phase transition point. In ID there 
is no superfluidity in the strict sense of the word, and no phase transition occurs, but 
the picture remains very similar to that of 2D and 3D case for a weakly interacting 
gas, as explained in the next subsection. A detailed discussion of the fluctuation 
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region, including accurate expressions for the critical points, will be presented later, 
in section [HI Meanwhile, here we want to utilize the fact that the results f|5.13f> - (|5.19[) 
allow one to make order-of-magnitude estimates of the non-perturbative fluctuation 
contributions to the thermodynamic functions. These estimates are quite important 
as they show the degree to which the perturbative results of the previous sections are 
inaccurate in the fluctuation region. As we will see, for some quantities the fluctuation 
corrections are smaller than the leading ideal-gas contributions. 

By continuity, the order-of-magnitude estimates for fluctuation contributions can 
be obtained from (|5.13[) - (|5.19[) by simply setting 77- ~ 1, while for the quantities 
themselves we can use their mean-field critical expressions. This way we arrive at the 
following results 



(A/Lt) fluc t (An) 



fluct 



7o /3 (d = 3), 
ln-^l/To) (d=2), (5.20) 
f (d=l). 



(Apjfluct (Ae)fl uct (As) fluct 
pes 

7o /3 (d = 3), 
7 ln(l/ 7o ) (d = 2), (5.21) 

7o /2 (d=l). 

For the superfluid and condensate densities we get an obvious answer that in both 
cases the fluctuation contributions are on the order of 100%. 



5.4- Specifics of ID system 

As long as a ID system is weakly interacting, i.e. jt -C 1, the notion of the 
order parameter field with well defined amplitude and the two-component (normal 
+ superfluid) description remain physically meaningful. Since superfluidity is a 
topological phenomenon, it can be destroyed only by topological defects — phase slips. 
At 7t < 1 the phase slips are rare events which do not contribute significantly to 
the local thermodynamic quantities; in correlation functions, phase slips show up only 
at length scales much larger than 1/fco, where their effect can be described at the 
hydrodynamic level. 

When the temperature reaches the characteristic scale of 

-70-^-70^, (5.22) 

the parameter 77- becomes of order unity. Apart from the lack of a genuine phase 
transition, the physics in the temperature range T ~ is close to that of the 

fluctuation region in 2D and 3D. Only the long- wave classical-field subsystem of the 
original quantum-field experiences strong non-linear fluctuations. This leads to non- 
perturbative contributions to the system thermodynamics which are universal for all 
weakly interacting one-dimensional U(l) systems, and can be described by universal 
scaling functions in direct analogy with the fluctuation regions in 2D and 3D p"5l ITT] . 
At T ^> Tg^t , the low-momenta part of the classical-field component gets depleted, 
so that the non-linearity of interactions becomes weak and accurately accounted for 
within the normal-gas mean-field picture. 
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6. Off-diagonal correlations 

Bcliacv's diagrammatic technique allows one to calculate correlation functions up to 
distances much larger than the correlation radius ro ~ 1/ko- However, addressing 
the asymptotic long-range behaviour of off-diagonal correlation functions requires 
special techniques properly accounting for long-wave fluctuations of the Goldstone 
mode, i.e. the phase of the superfluid order parameter. Since the diagrammatic and 
hydrodynamic description overlap, one can straightforwardly extend the diagrammatic 
calculation of off-diagonal correlators obtained up to large enough distances, r ^> tq to 
arbitrarily large r's. We proceed in the spirit of Popov's hydrodynamic approach |19j . 
but with a significant simplification that the hydrodynamic treatment takes place on 
top of Bcliacv's techniques, and thus does not require any additional modifications. 
The simplification is possible due to the fact that at large distances only phase 
fluctuations are important and their effect can be factored out as 

V> = ^c i<E \ (6.1) 

where, under coarse-graining up to momentum k\ , 

$(r,r) = ]T e ikr $(k,r) (6.2) 

k<ki 

is the long-wave part of the phase field. As before, the momentum cutoff k\ <C fco 
should be small enough to guarantee (i) the statistical independence of the fluctuations 
of 4> from the fluctuations of ip, and (ii) the Gaussian character of phase fluctuations 
in the hydrodynamic regime. (A particular value of ki is not important and does not 
appear in final expressions.) 

In view of (i) and (ii), we have for an m-particlc correlator, 

K, = ( V>1 V>2 ^3 i>4 ■ ■ ■ i>2m-l V>2m ) ) ( 6 - 3 ) 

where the subscripts label the space-time variables. Using (|6.2p it can be written as 
K = ICc- A , (6.4) 

with 

1 ^ m 2 

A = o( (Et- 1 )^) )< ( 6 - 5 ) 

J'=l 

where IC is obtained from IC by substituting ipj — > for all j's. 

At distances | — r s | 3> 1/fco the correlator /C remains constant, and the remaining 
dependence on distance is exclusively due to the fluctuations of 4>. [The correlator IC 
does depend on the coarse-graining momentum k\ and this dependence is crucial for 
compensating the fci-dependence of the function A; in fact its origin is directly related 
to A when coarse-graining is stopped at larger momenta and then taken further to k\ .] 
The independence of IC on coordinates and times in the asymptotic limit immediately 
leads to the relation 

IC(X) = K(X')e A ^- A & (6.6) 

between correlators at different sets of variables, X and X', provided both are in the 
asymptotic region. Now if X' is within reach of diagrammatic expansions, then (|6.6p 
allows one to "extrapolate" IC(X') to IC(X) at arbitrary X in the asymptotic domain. 
By the structure of (|6.6|) . the difference A(X') — A(X) is independent of hi, as opposed 
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to individual A's. Here we assumed for simplicity that all variables in X and X' arc 
large enough. The same analysis is readily adapted for cases when only some of the 
coordinates in the m-particle correlation function arc in the asymptotic regime; only 
these coordinates have to be mentioned in (|6.6p . 

The Gaussian character of the field $ implies, see (|6.5[) . that 



A(X) - A(X>) = ]T (-i)>+i E(v sj ,r sj ;v' sj ,r' 8j ) , (6.7) 

where r S j = r s — Tj, r s j = t s — Tj (the same for primed variables), and 
H(r, t; t' ,t') = ( $(r, r)$(0, 0) - $(r', t')$(0, 0) ) 
= T J2 [c i(kr ^ r) - e i(kr '-« r '>] ( |$(k,0| 2 ) • (6.8) 

We see that for the extrapolation of correlation functions to arbitrarily large distances 
one needs to know only E(r, r; r', r') which is defined through the average ( |$(k, £)| 2 ). 
The latter is readily found from Popov's hydrodynamic action 

S = |£[K/m)fc 2 + < 2 ]|$(k,£)| 2 , (6.9) 

leading to 

(|$(k,C)| 2 ) = [K/m)fc 2 + < 2 ]^ 1 . (6.10) 

In fact, non-zero frequencies (accounting for the quantization of the fluctuations 
of the phase) are relevant only at T <C nil , meaning that the parameter k in (|6.10D 
can be taken in its T = limit, 

AtT> nil, only the ^ = term should be left in (|6Tg|) . 
Let us consider now the single-particle density matrix 

p(r) = (^(r)^(O)). (6.12) 

Taking into account that p(r) = n + G(r = 0, r = —0) — G(r, t = —0), from the 
diagrammatic expressions for G we have (at small and intermediate r's) 

^) = »-£(l-c-)[i§^(l + 2» E )-i]. (6.13) 

We can use this expression in any dimension at distances significantly exceeding 1/fco- 
In 2D at finite temperature and in ID at any temperature, the expression (|6.13l) 
ultimately becomes inaccurate, and we have to rely on the above-described procedure 
of extrapolation. In the asymptotic regime we have 

p(r) = p(r')c- s( - r - r '^ , (6.14) 

where, at T -C nil, 

-( r > r ) = -l^[ e " e J ^ . ( 6 - 15 ) 

k 

with c = \J n s /mn the sound velocity at T — (here n s = n), while at T > nil : 

mT 



ml r =i._< :i— 1 1 
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Within the relevant orders, equation (|6.13[) corresponds to the expansion of the 
exponential in (|6.14[) in powers of H, up to the term oc S included. This immediately 
suggests that within the same accuracy one can exponentiate S-tcrms in (|6.13[) to 
extend its domain of applicability to much larger distances. The other advantage of 
proceeding this way is having a physical definition of the quasi-condensate density as 
the amplitude of the order parameter field in the long- wave limit [13j . The equations 
that follow have the same accuracy as (|6.13p though they do contain artificial higher- 
order terms which arise from factorizing the correlation function. More specifically, we 
single out terms in (|6.13p which on large distances reproduce S, and then exponentiate 
them: 

p(r) = [n - d(r) - nA(r)] — ■+ p(r) e" A M , (6.17) 
where p(r) = n — d(r), and 

d(r) =E( X - ^0 ^ [ E ~l + n + {E - ® Ne ] ' (6 ' 18) 

k 

Pr/i „ikr\ E - e 



A « - -~ E t 1 - eikr ) w^ 1 + 2Ne] ■ (fU9) 

k 

We have added and subtracted [1 + 2N E ] [pe/2E 2 ] to (|6.13[) to ensure that (i) both 
p(r) and A(r) are free from ultra-violet and infra-red divergencies in all dimensions, 
and (ii) that only small momenta k < ko contribute to the phase correlator A(r). 

By comparing (|6 . 19[) and (|6 . 1 5|) we see that they coincide at large distances up to 
leading terms. It means that in 2D and in ID at zero temperature, equation (|6 . 1 T[) can 
be trusted up to exponentially large scales, while in ID at finite temperature T > \p,\ 
it works at least up to distances ~ n s /mT. 

We are now in position to define the quasi-condensate density as the limiting 
value of p(r — > oo): 

-E- e + p, 



+ (E-p)N E . (6.20) 



The physical meaning of this relation is the amplitude squared of the order parameter 
field at large distances. There is a certain degree of freedom in attributing terms 
which do not result in the power-law or exponential decay of p(r) to d(r) or to A(r). 
The idea behind our choice is three-fold: (i) equations (|6. 171) . (|6.18p . and (|6.19[) are 
valid as written in all spatial dimensions; (ii) at large distances A(r) has the structure 
of hydrodynamic phase correlations; (iii) the final expressions have the simplest form 
possible within the same accuracy. It is also worth mentioning that n qc is a quantity 
which controls all physical processes happening in the WIGB at short distances at low 
temperature, e.g. m-body recombination rates. In all spatial dimensions it plays the 
same role as the condensate density in 3D system as long as one is interested in length 
scales not much larger than the healing length [6l [13] . 
Finally, the condensate density is defined from 

n = n qc e- A(oo) . (6.21) 

It is not accidental that the ultimate long-wave length property of the supcrfluid 
system is determined last. It was always an unpleasant feature of numerous mean- 
field treatments that the crucial parameter determining physics at short scales was 
linked to uq. 
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7. Normal region 

As was mentioned in the Introduction, we assume that the temperature is not very 
high, so that the condition (|1.2j) is preserved and the quantity U remains the only 
parameter characterizing the interaction between the particles. 



7.1. Thermodynamic functions 



In the normal region, where there are no anomalous correlators, we have to deal 
only with the Dyson equation for the Green's function G in terms of self-energy S n . 
Within the leading-order approximation, Six = 2n£7 (in 2D this formula implies an 
adequate choice of the parameter eo = £q(T) for the effective interaction U, discussed 
in subsection l7.2p . and the expression for G yields the self-consistent relation for the 
number density, 

-l 



i 



with 



~e(k) = e(k) - m = e(k) + 
For the pressure we find 



(7.1) 
(7.2) 



V 



» d» 
n — d/i = 

-oo d M 

nd/z + n 2 U 



n ( 1 + 2U— I d/i 

-oo V d/i 



Using 



we finally arrive at 



n U - T 



-m/T 



(7.3) 



(7.4) 



(7.5) 



Equations (|7.ip . (|4.3p . and (|7.5p specify n, /i, and p as functions of (T,p,). Utilizing 
(|4.29p . with x = p., and then (|4.30p . we find for the entropy 



and energy 



= J2 \ ^-N- e -ln [l- e -m)/T 



(7.6) 



(7.7) 
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Figure 11. Second-order contributions to the self-energy in the normal regime, 
including the sunrise diagram. 



7.2. Effective interaction in the normal regime in 2D 

In the quasi-condensate region, the parameter eo defining the effective interaction U 
can vary over a wide range of values thanks to the retained second-order correction 
that produces counter-terms compensating for the arbitrary choice. Since in the 
normal region we confine ourselves to the first-order expression for the self-energy, 
S = 2nU , we need to choose the value of eo which minimizes the omitted second-order 
contributions originating from the sunrise diagram, see figure QT] Due to momentum 
independence of the pseudo-potential line, the third (fourth) diagram is identical to 
the first (second) one; we thus consider only the first two diagrams and multiply the 
result by a factor of 2. The parameter eo can be much smaller or much larger than 
T . In either case the leading term comes from the logarithmic ultraviolet contribution 
from the product of two propagators, yielding the result 



Comparing the first two terms in the brackets with (|3.51|) and (|3.52|) . we see that 
if we take the value of eo substantially away from T, then the leading second-order 
correction to the expression S = 2nU amounts to renormalizing the value of the 
pseudo-potential U in such a way that it corresponds to eo ~ T. This brings us to 
the conclusion that eo ~ T is the optimal choice for eo, in which case omitting the 
diagrams shown in figure [TT1 is justified by the parameter mil <C 1. 

7.3. Expansion parameter 

At temperatures above the fluctuation region where the rcnormalizcd chemical 
potential remains small \p,\ <C T the expansion parameter is defined by the infra- 
red behaviour of the zero-frequency lines. The dimcnsionless factor, jt, associated 
with adding an extra full interaction vertex to a diagram is given by (|5.8p . but now 
with G ~ T/e(fc), resulting in the estimate 




(7.8) 



k>k T 




(7.9) 




(7.10) 
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where /2 fluc is the size of the fluctuation region in terms of the chemical potential. 
At temperatures much higher than the degeneracy temperature T c ~ n 2 / d /m where 
jj, m — Tln(nAy), all corrections are small in the parameter 

(7 ' n) 



8. Pseudo-Hamiltonian 



In Section [4] we have obtained the thermodynamics of the system in the quasi- 
condensate region, see (|4.17|) - (|4.18|) , (|4.28[> . (|4.31|) , and (|4.32|) which specify n, /i, 
p, s, and e as functions of two independent variables, fl and T, in the form of integrals 
involving the /i-dependent quasiparticlc spectrum, E(k), given by (|4.6|) . So far, we did 
not discuss the physical meaning of the function E(k), since this was not necessary for 
our derivations. On the other hand, equation (|4.3ip for the entropy clearly suggests 
that, within our approximation, the system is equivalent to a system of non-interacting 
bosonic quasiparticles with dispersion E(k), described by the Hamiltonian 

H = E + Y, E(k)n k , (8.1) 

k 

where is the quasiparticlc occupation number operator and E$ is some constant. 
The natural question is then whether the Hamiltonian (|8.ip can be used for calculating 
other thermodynamic quantities within the non-interacting quasiparticlc gas picture. 
By direct comparison with (|4TT7|) - (|4TT5|) . (|4T2"8"]1 . P~3"T|) . and (|432"|) , it can be shown 
that this is indeed the case and one can write 

E- f iN= (H) =E +Y, E{k)N E , (8.2) 

k 

Q = -pV = -Tin (Trc~ H/T ) . (8.3) 
It is important, however, to remember that Eq is a function of both fl and T, satisfying 

E = -^j + 2^n' -Un' 2 + l[ d) (ra). (8.4) 

In view of this temperature dependence, we refer to (|8.1[) as a pseurfo-Hamiltonian, 
rather than an effective Hamiltonian. 

The pseudo-Hamiltonian description can be applied to the normal region as well. 
In this case, 

E = -n 2 U, (8.5) 

were n = n(p,, T) is specified by (|7.1[) . The /i-dependent spectrum of quasi-particles is 
now given by (|7.2[) . 



9. Fluctuation region 

In the fluctuation region the long-wave part of the classical-field component of the 
quantum field experiences strong (non-perturbativc) fluctuations. It is exclusively 
due to these fluctuations that the diagrammatic technique for the quantum field loses 
an expansion parameter. The special role of the classical-field component is evident 
from the fact that on approach to the fluctuation region the leading contribution to the 
expansion parameter 77- is associated with zero- frequency propagators, see subsections 



Beliaev technique for a weakly interacting Bose gas 



37 



15.21 and 17.31 The diagrammatic technique built on zero-frequency propagators (with 
perturbative parts of self-energies included) and the interaction represented by pseudo- 
potential U corresponds to the diagrammatic expansion of the Gibbs distribution 



Z = J C - H ^ T V^ (9.1) 
for the in momentum space truncated classical field 

^(r)= ^e ik ' r , (9.2) 



k<k> 

with the Hamiltonian functional 




^|V<A| 2 + ||V'I 4 -mW 

2m 2 



dr. (9.3) 



Here k! <C kx = ymT is a truncation momentum (to be discussed later), and 
fi' = /i'(fc') is the reduced chemical potential obtained from the bare one by subtracting 
relevant perturbative contributions from all the harmonics with k > k', including the 
quantum ones. The bottom line is that describing the quantum gas in the fluctuation 
region reduces to solving the classical- field problem (|9.ip - (|9.3p . in its fluctuation 
region. 

In al = 2,3, numerical solutions to the problem (|9.1[) - (|9.3p — in terms of scaling 
functions for thermodynamic quantities — are available in the literature and have 
already been used to accurately describe weakly interacting 2D and 3D quantum 
gases in the fluctuation region [T51 [lj)J HZ] ■ The same is possible in ID, but to the best 
of our knowledge it has not been done so far. 

In the fluctuation region, an important quantity is the momentum k <C kx 
separating weakly and strongly coupled classical modes. For a quantitatively accurate 
description of the fluctuating classical-field sub-system, the momentum k' separating 
classical modes of interest from the rest of the modes has to be much larger than k. 
However, as long as we are interested in the order-of-magnitude estimates, it is safe 
and convenient to set k' ~ fc, so that all the modes we are dealing with in (|9.1j) - (|9.3[) 
are strongly coupled. To estimate k we note that with k' ~ k all three terms in the 
Hamiltonian (|9.3[) have to be of the same order of magnitude, that is 

k 2 /m ~ \p.\ ~ nil , (9.4) 

where 

n~£j^| a = X>' (9-5) 

k<k k<k 

is the long- wavelength contribution to the total density, and fi = fJ-'(k' — > k). For h we 
have h ~ k A n^ and since k is separating strongly coupled long- wave harmonics from 
slightly perturbed short-wave ones, the order-of-magnitude estimate for nr follows by 
continuity from the ideal system formula: 

T T 

nr. - —7 . (9.6) 

Substituting this back into fO )) -(|9.5 p yields 

k = (m 2 TU)^ , (9.7) 
h - (m d T 2 U d ~ 2 )^ , (9.8) 
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\fi\ ~ (m d T 2 U 2 )^ . (9.9) 

The estimates (|9.8|) - (|9.9[) imply the following parameterization of the grand-canonical 
equation of state in the fluctuation region and its vicinity: 

n = nW(T) + (m d T 2 U d - 2 )^ X id Hx) , (9.10) 
-- ^ d) (T) (gll) 



(m d T 2 U 2 )^ 

Here \W(x) is a dimcnsionlcss scaling function of a dimensionless scaling variable x. 
In 2D and 3D, the quantities ni d '(T) and ni d \T) are the critical values of density 
and chemical potential at a given temperature. In ID, where the phase transition 
is absent, one can set, without loss of generality, /i c = 0, and, correspondingly, 
Hc(T) = n(A* = 0,T). 

Similarly, supcrfluid and condensate densities — in the dimensions in which they 
are meaningful — are parameterized as 

n s =(m d T 2 U d - 2 )^fW(x) (d = 2,3), (9.12) 

and 

n = m 3 T 2 Uf (x) {d = 3) . (9.13) 

The functions X(x), f s {x), and fo(x) are universal for all weakly interacting U(l)- 
symmetric systems in the given dimension (no matter quantum or classical, continuous 
space or lattice); the numerical data for them is available for d = 2,3 in (T5J |T7] . By 
numerically solving the problem (|9.1|) - (|9.3p . it was established that, up to higher-order 
corrections in the parameter 70, [151 1161 [T7] 

nf D) = ni 0) (T) - Cm 3 T 2 U , C = 0.0142(4) , (9.14) 

nf D) = ^La(X). e = 380±3, (9.15) 
In \mu J 

with ni°^ (T) the critical density for Bose-Einstein condensation in an ideal 3D gas. 
Analogous relations for /j,c (T) read [THl [THl E] 

^ = 2Un?HT) + ^U°4?m). (9.16) 

7T 2 VVm 3 riW 

(2 D) = rnTU ^ /_^\ = 13 2 ± Q 4 

Note that the leading terms in the above relations, as well as ordcr-of-magnitudc 
estimates for the values of sub-leading terms, are readily obtained from the condition 
7t ~ 1; the accurate numerical treatment of the problem ()9.1[) - (|9.3[1 is necessary to 
fix the values of the dimensionless constants. 

By re- writing (|9.15| . (|9.14|) in terms of the critical temperature as a function of 
density, we get 

rr } = ; m , ^ 

m m{t, / mu ) 

=T^(1 + C an 1 / 3 ), C = 1.29 ± 0.05 , (9.19) 
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Figure 12. (Colour online) Density matrix for the two-dimensional Bosc- 
Hubbard model (fTo7T] l at U/t = 0.25, T/t = 0.5, and n = 0.5 with periodic 
boundary conditions on a square lattice lattice with 200 X 200 sites. Theoretical 
lines were obtained from J 6 . 13 ft (dashed red line) and 116.171 (solid black line). 



where is the critical temperature of the ideal 3D gas. Now it is easy to estimate 
k. The relevant temperature for having strong fluctuations is given by 

T fl uct«T c (d = 2,3), T nuc t ~ nUho (d = 1) . (9.20) 

The estimate (|9.7p at T = Xh uct then yields 

(d = 3) , 



2/3 

7o 

k/k Tnuct ~ { 7o 1/2 (d = 2) , 
7o /2 (d=l). 



(9.21) 



We see that fcr fluct in all cases, which justifies the statement that the physics of 
the fluctuation region is that of a (strongly interacting) classical field. 

Putting aside n s and no, which are most sensitive to non-perturbative fluctuations 
of the classical field ip, the next two quantities which are sensitive to fluctuations, 
especially in ID and 2D, are the density and the chemical potential — see estimates 
(|5.20[) . For other quantities — as is seen from (|5.2ip — neglecting fluctuation corrections 
will not result in a significant error. It is also important that up to a few dimensionlcss 
constants characterizing sub- leading contributions to critical values of p, e, and s, the 
fluctuation corrections to these quantities are expressed in terms of the same function 

(x) and its integral. Indeed, equation (|4.19[) implies the following relation for p: 



P 



: pW(T) + nW(T)[ M . 



/'c 



+ (m 2d T i U d )^ / A (d V)dx' 



(9.22) 



and then with (|4.29[) - (|4.30[) one obtains similar relations for s and e. 



10. Comparison with numerical results 



To see how accurate our description is for the density matrix we performed Monte 
carlo simulations of weakly interacting two- and onc-dimensional systems where phase 
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Figure 13. (Colour online) Density matrix for the one-dimensional Bose gas with 
na = 8 at temperature mT/(2irn 2 ) = 0.001 (black circles) and mT/ (2nn 2 ) = 
0.004 (red squares) on a circle of length L/a = 25. Theoretical curves are 
represented by 116.171 1. 1 16. 181 1. and 1 16. 191 1 with the dashed black line for the lower 
temperature and the solid red line for higher one. 

fluctuations are the strongest. We deliberately aim at systems with substantial decay 
of the density matrix to see the difference between the perturbative treatment and 
exact solution. In figure rT2] we present data for the two dimensional square-lattice 
system described by the Bosc-Hubbard Hamiltonian 

H = -tJ2 (VrW-'r + h - C ) + Y - 1] - (10-1) 

<rr'> r 

where t is the hopping amplitude between the nearest neighbor sites. Model (| 1 . 1[) 
was simulated for U/t = 0.25, T/t = 0.5 and filling factor n = 0.5 on a lattice with 
L x L = 200 x 200 points by using the Worm Algorithm approach [31] . Formally, all 
expressions in this manuscript relating energy spectrum, particle density, and density 
matrix remain valid for an arbitrary dispersion relation with e(k) — ► k 2 /2m at small 
momenta. Moreover, the theory is supposed to work as written for systems with 
periodic boundary conditions provided the sums over momenta are understood as 
L~ d X)k^o> where k = lnxv/L and n is a <i-dimensional integer. To suppress statistical 
noise, simulation data for p(r) were collected to spherically symmetric bins. We are 
using exactly the same procedure to present theoretical data. Thus, the only difference 
between the essentially exact (up to statistical errors) simulation data and the theory 
is due to the finite value of U. Even for relatively large values of U and temperature 
T ~ 0.6T C we observe a remarkable accuracy of (|6.13|) . One can get an idea of the 
systematic theoretical error by comparing curves derived from (|6.13[) and (|6.17[) . [The 
latter is obtained by exponentiating the phase correlator at large distances]. 

In figure [13] we apply our theory to the Lieb-Liniger model of the one-dimensional 
Bose gas with contact interactions. Since exact correlation functions at finite 
temperature are not known, we performed Monte carlo simulations using recently 
developed techniques for continuous space systems [321 133] . It is convenient to 
introduce the on-dimensional characteristic length Iq = 2 /mil as the unit of length. 
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Figure 14. Energy per particle in units of T c for the three-dimensional Bose 
gas with na 3 = 10 — 6 vs. temperature, in the quasi-condensate region. Squares 
are quantum Monte Carlo results, solid line is the theoretical curve 1 14.321 . 

The limit of weak interactions is obtained then by considering large densities na 1, 
and we have chosen nlo = 8 corresponding to the power law-decay of off-diagonal 
correlations at zero temperature with the exponent a = l/w-\/2nlo = 1/Att. The 
upper curve in figure [T51 shows results for low temperature mT/(2nn 2 ) = 0.001, or 
T/nU ~ 0.025 when thermal phase fluctuations are barely influencing the data for the 
simulated system size L/Iq — 25. When temperature is increased to mT/(2im 2 ) = 
0.004 (lower curve in figure [T5]) we clearly see the bimodal decay of the density matrix 
and the crossover from the power-law to exponential decay. Contrary to the 2D 
case, equation (|6 . 1 T[) captures the actual behaviour better than (|6 . 13[) which is not so 
surprising given the large suppression of p(r) at large distances. 

Finally, figures [17] and [18] show the agreement of Monte Carlo results for the 
energy (squares) and the pressure (triangles) at temperatures T > T c with theoretical 
curves from (|7.7|) and (|7.5|) (solid lines). 

In the following we compare thermodynamic functions obtained in Sections [4] and 
[7] with exact quantum Monte Carlo results of a weakly interacting three-dimensional 
system. Let us start with the case of a homogeneous system with a small parameter 
na? = 10~ 6 . In figures [T4l and [TBI we compare Monte Carlo results for the energy 
(squares) per particle and pressure (triangles) at temperatures T < T c , with the 
theoretical curves from (|4.32[) and (|4.28[) (solid lines). The analytical expressions are 
in very good agreement with numerical results up to temperatures T ~ T c and there 
is no need to switch to the description in terms of universal functions |15| . This is 
consistent with the estimate (|5 . 2 1 1) showing that fluctuation contributions to energy 
and pressure are very small. The same does not apply to the superfluid and condensate 
densities for which the perturabtive expansion in U is not valid as fluctuations of the 
order parameter cannot be neglected on approach to the critical temperature. In 
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Figure 15. Pressure in units of l/2ma 2 for the three-dimensional Bose gas 
with no 3 = 10 — 6 vs. temperature in the quasi-condensate region. Triangles are 
quantum Monte Carlo results, solid line is the theoretical curve 114.281 1. 




Figure 16. Quantum Monte Carlo results for superfluid fraction (squares) 
and condensate fraction (triangles) for the three-dimensional Bose gas with 
no 3 = 10 -6 . Solid and dashed lines are a combination of theoretical expressions 
from 114.591 1 for the superfluid density and 1 14. 71 for condensate density with 
preexisting classical Monte Carlo results 1151 for the fluctuation region (the arrow 
is pointing to the gap in the curves where the two results are matched). 
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Figure 17. Deviation of the energy per particle (in units of T c ) from the classical 
law for the three-dimensional Bose gas with na 3 = 10 — 6 vs. temperature, in 
the normal region. Squares are quantum Monte Carlo results, solid line is the 
theoretical curve i7.7l . 




Figure 18. Deviation of the pressure (in units of l/2ma 2 ) from the classical 
law for the three-dimensional Bose gas with na 3 = 10 — 6 vs. temperature, in 
the normal region. Triangles are quantum Monte Carlo results, solid line is the 
theoretical curve 17.51. 
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figure (fTB|) squares and triangles are the results of quantum Monte Carlo simulations 
for the supcrfluid and condensate fraction, respectively. Solid and dashed lines are 
the theoretical expressions from (|4.59p for the supcrfluid density and (|4.7[) for the 
condensate density, respectively, combined with preexisting classical Monte Carlo 
results [15] for the fluctuation region. At a temperature T > 0.6T C we switch to 

the universal description, as 7^ u1 ^ [see (|5.12|) ] becomes of the order of unity. For 
T ~ T c the remaining discrepancy between classical and quantum Monte Carlo results 
is due to finite interaction strength. One should not be misled by the small value of 
the gas parameter in the fluctuation region since the proper parameter controlling the 
size of the fluctuation region in temperature, AT/T C , and the accuracy of the universal 
description involves a large numerical prefactor and is rather na 3 /10 -5 , see |15j . 

11. Conclusions 

We have shown that Bcliaev's diagrammatic technique regularized by adding a small 
term explicitly breaking the U(l) symmetry to the Hamiltonian yields a simple and 
controllable description of the weakly interacting Bose gas in any dimension. This 
approach is especially convenient for obtaining thermodynamic functions that are 
not sensitive (up to higher-order corrections) to the long-range phase fluctuations. 
The symmetry breaking term introduces a small gap in the otherwise gapless 
Goldstone mode and suppresses long-range fluctuations of the phase of the order 
parameter, thereby introducing a genuine condensate density and thus substantially 
simplifying the description of lower-dimensional systems. In this sense, the symmetry- 
breaking trick is an alternative to Berezinskii's finite-size trick and Popov's special 
(hydrodynamic) treatment of long- wave parts of the fields. At the end of the 
calculation, the symmetry breaking term is set to zero. 

Another important feature of our approach is that it completely avoids such 
notions as seed- or quasi-condensate for its construction. The quasi-condensate density 
defined as the modulus squared of the order parameter field at distances larger then 
the healing length is calculated as the long-wave property of the theory, rather than 
serving as an input parameter. It is only natural to have the theory which handles 
the short-range physics first and uses it for dealing with the low-energy physics next. 
As far as long-range off-diagonal correlations are concerned, the approach produces 
accurate answers for the off-diagonal correlation functions up to distances where 
further evolution of the correlators is controlled by generic hydrodynamic relations, 
and thus can be accurately extrapolated to arbitrarily large distances. 

We confined our analysis to the most typical and non-trivial (in 2D and 3D) 
case of dilute gas, when the size of the potential Rq is much smaller than the distance 
between the particles. In the (Boltzmann) high-temperature regime, we have restricted 
ourselves to low enough temperatures at which dc Broglic wavelength remains much 
larger than R$ and the answers do not depend on the details of the pair potential. 
Under the two above-mentioned conditions, the effective inter-particle interaction is 
described by a single parameter U. Technically, generalization of the theory to the case 
of Rq > n~ x / d is straightforward. In this case, the weakness of interaction literally 
means Born interaction potential, and complete summation of ladder diagrams in 2D 
and 3D becomes irrelevant. The description of a normal gas at \t < n~ x / d is a 
well-studied topic that goes beyond the scope of the present paper. 

The pair of variables fl = jjL — 2nU and T form the most convenient set 
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of parameters determining the rest of thermodynamic quantities, including the 
condensate density, if any. The structure of the answers corresponds to a picture 
of non- interacting bosonic quasiparticles with the Bogoliubov spectrum (|4.6[) . on top 
of the 'vacuum' with temperature-dependent energy. This allows one to introduce the 
bilinear bosonic pscudo-Hamiltonian (|8.ip with Bogoliubov spectrum. We argue that 
within the bilinear pseudo-Hamiltonian ansatz our results cannot be further improved 
at least at T > nil . In the T <C nil limit, when thermodynamics is exhausted by dilute 
non-interacting phonons, using the improved value of the sound velocity obtained 
directly from the zero-temperature compressibility, yields more accurate answers for 
thermodynamic quantities, but is inconsistent with retaining the Bogoliubov form of 
the spectrum for all the momenta. We also show that any attempt to improve the 
self-consistent mean-field description of the superfluid region, aimed at getting rid of 
the spurious first-order phase transition, is senseless as long as the fluctuation-induced 
non-perturbative shift of the critical temperature is not accounted for. In fact, there 
is no need to extrapolate the pscudo-Hamiltonian description to the fluctuation region 
around the critical temperature in view of the availability of the numerical answers 
for universal scaling functions describing all the weakly interacting U(l) systems at 
\T — T c \/T c <C 1. With these data, the two analytical descriptions — in the normal and 
superfluid phases — are readily connected across the fluctuation region. By comparing 
our results to first-principles Mote Carlo data, we observe that the most vulnerable 
region of parameters is \T — T c \/T c -C 1, where the description is accurate only when 
the gas parameter an 1 / 3 is as small as ~ 10~ 2 . 
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